In [1]:
%pylab inline
Ecuación a resolver $$ \frac{\partial u}{\partial t} = \nu\frac{\partial^2 u}{\partial x^2} $$
Discretizando
$$ \frac{\partial u}{\partial t} \approx \frac{u_{i}^{j+1} - u_{i}^{j}}{\Delta t} $$$$ \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1}^j - 2 u_{i}^{j} + u_{i-1}^{j}}{(\Delta x)^2} $$En el esquema de diferencias finitas $$ \frac{u_{i}^{j+1} - u_{i}^{j}}{\Delta t} = \nu \frac{u_{i+1}^j - 2 u_{i}^{j} + u_{i-1}^{j}}{(\Delta x)^2} $$
$$ u_{i}^{j+1} = \nu\alpha u_{i+1}^{j} + (1 - 2\nu\alpha)u_{i}^{j} + \nu\alpha u_{i-1}^{j} $$donde $\alpha = \Delta t/ (\Delta x)^2$.
In [25]:
n_x = 100
n_t = 100
nu = 0.3
sigma = 0.2 #parámetro que asegura que \alpha\nu < 0.5
x = linspace(0, 1.0, n_x)
dx = x[1]-x[0]
dt = sigma*dx**2/nu #dt se define utilizando sigma
alpha = dt/dx**2
print dt
u = zeros(n_x)
u[where((x<0.66) & (x>0.33))] = 1.0
In [26]:
plt.plot(x,u)
plt.title('Condicion inicial')
plt.xlabel('x')
plt.ylabel('u')
Out[26]:
In [27]:
for n in range(n_t): # loop over time
u_past = u.copy()
for i in range(1,n_x-1): #loop over space
u[i] = nu * alpha * u_past[i+1] + (1.0 - 2.0*nu*alpha)*u_past[i] + nu*alpha*u_past[i-1]
In [28]:
plt.plot(x,u)
plt.title('Condicion despes de $\Delta t=100$')
plt.xlabel('x')
plt.ylabel('u')
Out[28]:
In [ ]: